
#include <algorithm>
#include "deformationMesh_rbfguided.h"
#include "weightsScheme.h"

//~ namespace MLS {

DeformationMesh_RBFguided::DeformationMesh_RBFguided() {
   DeformationMesh::DeformationMesh();
   
   rbf = 0;
}

DeformationMesh_RBFguided::DeformationMesh_RBFguided(Mesh *m) {
   
   DeformationMesh::DeformationMesh();
   assert(m);
   model = m;
   rbf = 0;
}

DeformationMesh_RBFguided::~DeformationMesh_RBFguided() {
   
}
   
void DeformationMesh_RBFguided::setting(Constraints *constraints) {
   
   assert(model);
   //~ DeformationMesh::setMesh(model); ???????????????
   DeformationMesh::setting(constraints);
   
   //~ // Compute the weights table
   //~ if (WeightsTable != 0) delete []WeightsTable;
   
   //~ unsigned int nminusone = std::count(constraints->ids.begin(), constraints->ids.end(), -1);
   //~ if (nminusone == 0) // If all constraints are associated with joints
      //~ WeightsTable = MLS::compute_WeightsTable(model, constraints);
   //~ else if (nminusone < constraints->ids.size()) // If constraints are represented a trihedron
      //~ WeightsTable = MLS::compute_WeightsTable_Trihedrons(model, constraints);
   //~ else
      //~ WeightsTable = 0;
   
   //~ Deformation::initialize();
   
   cout<<"RBF computing begin ..."<<endl;
   rbf = new RBFfunction;
   rbf->setting(constraints);
   cout<<"RBF computing begin ... done"<<endl;
}

void DeformationMesh_RBFguided::execute_deformation(Constraints *constraints) {
   
   rbf->compute_lambdas(constraints);
   d_constraintsPos = &constraints->positions; // positions of deformed constraints 
   Vertex_iterator vi = model->vertices_begin(); unsigned int i=0, iv=0;
   while (vi != model->vertices_end()) {
      
      Point3 qq = rbf->evalf(pstarList[i]);
      qstar[0] = qq[0]; qstar[1] = qq[1]; qstar[2] = qq[2];
      computeCorrelationMatrix_points(i); 
      eigenvalue = computeLargestRootOfQuarticEquation() + traceM; 
      
      compute_u();
      
      /*
      // Rotate the vertex position
      computeRotation(&vhatList[i][0]);
      float *p = &model->vpositions[vi->id()];
      p[0] = qstar[0]+rv[0];
      p[1] = qstar[1]+rv[1];
      p[2] = qstar[2]+rv[2];
      
      // Rotate the normal vector
      computeRotation(&model->i_normals[iv]);
      model->normals[iv  ] = rv[0];
      model->normals[iv+1] = rv[1];
      model->normals[iv+2] = rv[2];
      */
      
      ++vi; ++i; iv+=3;
   }
}

//~ } // end MLS namespace
